The question is how much G x E we see in the RNAseq data.

To get started on this I ran the file “ProcessCount.R”. For this the data was TMM normalized, voom transformed, and then a gtXenvironment model was fit in limma.

library(limma)
library(tidyverse)
── Attaching packages ─────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 2.2.1     ✔ purrr   0.2.4
✔ tibble  1.3.4     ✔ dplyr   0.7.4
✔ tidyr   0.7.2     ✔ stringr 1.2.0
✔ readr   1.1.1     ✔ forcats 0.2.0
── Conflicts ────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
library(stringr)
library(modelr)
library(broom)

Attaching package: ‘broom’

The following object is masked from ‘package:modelr’:

    bootstrap
library(multidplyr)

Get the limma objects

load("voom.int.Rdata")

How many genes show GxE?

generate vector of interaction coefficients and then topTable

colnames(fit.int$coefficients)
int.coef.n <- colnames(fit.int$coefficients) %>% str_which(":") 
topTable(fit.int,coef = int.coef.n)

Not working

let’s just do an ANOVA for each gene

counts <- as.tibble(counts.voom.gr$E) %>% 
  mutate(geneIndex=1:n())
counts[1:10,1:10]
counts.long <- counts %>% gather(key="name",value="expression", -geneIndex) %>%
  mutate(gt=factor(str_extract(name,"RIL_[0-9]*")), 
         env=factor(str_extract(name,"CR|UN")), 
         rep = factor(str_extract(name,"Rep[0-9]"))) %>%
  select(-name) %>%
  group_by(geneIndex) %>% nest()
counts.long
counts.long$data[[1]]

Now fit an anova and do a VCA for each gene in parallel

fitAnova <- function(df) aov(expression ~ gt*env,data=df)
getPvals <- function(fit) {
  pvals <- summary(fit)[[1]]$`Pr(>F)`
  names(pvals) <- dimnames(summary(fit)[[1]])[[1]] %>% str_trim()
  data.frame(t(pvals))
}
fitVCA <- function(df) {
  tmpvca <- anovaVCA(expression ~ gt*env,Data=as.data.frame(df))
  tmpvca$aov.tab[,"%Total"]
  }
  
cl <- create_cluster(cores=4)
Initialising 4 core cluster.
counts.long.test <- partition(counts.long[1:100,],geneIndex,cluster = cl)
group_indices_.grouped_df ignores extra arguments
cluster_library(cl,c("purrr","broom","VCA","stringr"))
cluster_copy(cl,fitAnova)
cluster_copy(cl,getPvals)
cluster_copy(cl,fitVCA)
system.time(
  counts.long.test <- counts.long.test %>% mutate(aov = map(data,fitAnova), 
                                                  glance = map(aov,glance),
                                                  pvals = map(aov,getPvals),
                                                  varcomp = map(data,fitVCA))
)
   user  system elapsed 
  0.039   0.041  17.603 
LS0tCnRpdGxlOiAiRyB4IEUgZXhwcmVzc2lvbiIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKVGhlIHF1ZXN0aW9uIGlzIGhvdyBtdWNoIEcgeCBFIHdlIHNlZSBpbiB0aGUgUk5Bc2VxIGRhdGEuICAKClRvIGdldCBzdGFydGVkIG9uIHRoaXMgSSByYW4gdGhlIGZpbGUgIlByb2Nlc3NDb3VudC5SIi4gIEZvciB0aGlzIHRoZSBkYXRhIHdhcyBUTU0gbm9ybWFsaXplZCwgdm9vbSB0cmFuc2Zvcm1lZCwgYW5kIHRoZW4gYSBndFhlbnZpcm9ubWVudCBtb2RlbCB3YXMgZml0IGluIGxpbW1hLiAgCgpgYGB7cn0KbGlicmFyeShsaW1tYSkKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoc3RyaW5ncikKbGlicmFyeShtb2RlbHIpCmxpYnJhcnkoYnJvb20pCmxpYnJhcnkobXVsdGlkcGx5cikKbGlicmFyeShWQ0EpCmBgYAoKR2V0IHRoZSBsaW1tYSBvYmplY3RzCmBgYHtyfQpsb2FkKCJ2b29tLmludC5SZGF0YSIpCmBgYAoKSG93IG1hbnkgZ2VuZXMgc2hvdyBHeEU/CgpnZW5lcmF0ZSB2ZWN0b3Igb2YgaW50ZXJhY3Rpb24gY29lZmZpY2llbnRzIGFuZCB0aGVuIHRvcFRhYmxlCmBgYHtyLCBldmFsPUZBTFNFfQpjb2xuYW1lcyhmaXQuaW50JGNvZWZmaWNpZW50cykKaW50LmNvZWYubiA8LSBjb2xuYW1lcyhmaXQuaW50JGNvZWZmaWNpZW50cykgJT4lIHN0cl93aGljaCgiOiIpIAp0b3BUYWJsZShmaXQuaW50LGNvZWYgPSBpbnQuY29lZi5uKQpgYGAKCk5vdCB3b3JraW5nCgpsZXQncyBqdXN0IGRvIGFuIEFOT1ZBIGZvciBlYWNoIGdlbmUKCmBgYHtyfQpjb3VudHMgPC0gYXMudGliYmxlKGNvdW50cy52b29tLmdyJEUpICU+JSAKICBtdXRhdGUoZ2VuZUluZGV4PTE6bigpKQpjb3VudHNbMToxMCwxOjEwXQpgYGAKCmBgYHtyfQpjb3VudHMubG9uZyA8LSBjb3VudHMgJT4lIGdhdGhlcihrZXk9Im5hbWUiLHZhbHVlPSJleHByZXNzaW9uIiwgLWdlbmVJbmRleCkgJT4lCiAgbXV0YXRlKGd0PWZhY3RvcihzdHJfZXh0cmFjdChuYW1lLCJSSUxfWzAtOV0qIikpLCAKICAgICAgICAgZW52PWZhY3RvcihzdHJfZXh0cmFjdChuYW1lLCJDUnxVTiIpKSwgCiAgICAgICAgIHJlcCA9IGZhY3RvcihzdHJfZXh0cmFjdChuYW1lLCJSZXBbMC05XSIpKSkgJT4lCiAgc2VsZWN0KC1uYW1lKSAlPiUKICBncm91cF9ieShnZW5lSW5kZXgpICU+JSBuZXN0KCkKY291bnRzLmxvbmcKY291bnRzLmxvbmckZGF0YVtbMV1dCmBgYAoKTm93IGZpdCBhbiBhbm92YSBhbmQgZG8gYSBWQ0EgZm9yIGVhY2ggZ2VuZSBpbiBwYXJhbGxlbApgYGB7cn0KZml0QW5vdmEgPC0gZnVuY3Rpb24oZGYpIGFvdihleHByZXNzaW9uIH4gZ3QqZW52LGRhdGE9ZGYpCgpnZXRQdmFscyA8LSBmdW5jdGlvbihmaXQpIHsKICBwdmFscyA8LSBzdW1tYXJ5KGZpdClbWzFdXSRgUHIoPkYpYAogIG5hbWVzKHB2YWxzKSA8LSBkaW1uYW1lcyhzdW1tYXJ5KGZpdClbWzFdXSlbWzFdXSAlPiUgc3RyX3RyaW0oKQogIGRhdGEuZnJhbWUodChwdmFscykpCn0KCmZpdFZDQSA8LSBmdW5jdGlvbihkZikgewogIHRtcHZjYSA8LSBhbm92YVZDQShleHByZXNzaW9uIH4gZ3QqZW52LERhdGE9YXMuZGF0YS5mcmFtZShkZikpCiAgdG1wdmNhJGFvdi50YWJbLCIlVG90YWwiXQp9CgpjbCA8LSBjcmVhdGVfY2x1c3Rlcihjb3Jlcz00KQpjb3VudHMubG9uZyA8LSBwYXJ0aXRpb24oY291bnRzLmxvbmcsZ2VuZUluZGV4LGNsdXN0ZXIgPSBjbCkKY2x1c3Rlcl9saWJyYXJ5KGNsLGMoInB1cnJyIiwiYnJvb20iLCJWQ0EiLCJzdHJpbmdyIikpCmNsdXN0ZXJfY29weShjbCxmaXRBbm92YSkKY2x1c3Rlcl9jb3B5KGNsLGdldFB2YWxzKQpjbHVzdGVyX2NvcHkoY2wsZml0VkNBKQoKc3lzdGVtLnRpbWUoCiAgY291bnRzLmxvbmcgPC0gY291bnRzLmxvbmcgJT4lIG11dGF0ZShhb3YgPSBtYXAoZGF0YSxmaXRBbm92YSksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZ2xhbmNlID0gbWFwKGFvdixnbGFuY2UpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcHZhbHMgPSBtYXAoYW92LGdldFB2YWxzKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHZhcmNvbXAgPSBtYXAoZGF0YSxmaXRWQ0EpKSAlPiUKICAgIHNlbGVjdCgtZGF0YSkKKQoKc2F2ZShjb3VudC5sb25nLGZpbGU9Ikd4RV9tb2RlbHMuUmRhdGEiKQpgYGAKCmBgYHtyfQoKYGBgCgoK